

import gdal

def checkRasterSupport(layer1,layer2):
    """ test if the support of two rasters is good enough
    to do map algebra.

    raster must have same number of rows and columns

    difference in xmax < half cell width
    difference in xmin < half cell width
    difference in ymax < half cell height
    difference in ymin < half cell height

     the half cell size tolerance allows for floating point precision and a lot more...

    return (True,extent) or (False,0)
    """

    # return this if it fails
    notGood = (False,0)

    # check columns:

    if layer1.width() != layer2.width():
        return notGood
    if layer1.height() != layer2.height():
        return notGood

    # compute cell widths
    e1=layer1.extent()
    w1=e1.width() / layer1.width()
    h1=e1.height() / layer1.height()
    
    e2=layer2.extent()
    w2=e2.width() / layer2.width()
    h2=e2.height() / layer2.height()

    w=min(w1,w2)
    h=min(h1,h2)
    
    if abs(e1.xMax()-e2.xMax()) > w:
        return notGood
    if abs(e1.xMin()-e2.xMin()) > w:
        return notGood
    if abs(e1.yMax()-e2.yMax()) > h:
        return notGood
    if abs(e1.yMin()-e2.yMin()) > h:
        return notGood

    # compute average extent and return, in case the two are slightly different:
    xMax=(e1.xMax()+e2.xMax())/2
    yMax=(e1.yMax()+e2.yMax())/2
    xMin=(e1.xMin()+e2.xMin())/2
    yMin=(e1.yMin()+e2.yMin())/2
    
    return (True, (xMin,yMin,xMax,yMax))

    

def layerAsArray(layer):
    """ read the data from a single-band layer into a numpy array.
    Only works for gdal layers!
    """
    
    gd = gdal.Open(str(layer.source()))
    array = gd.ReadAsArray()
    return array

def writeGeoTiff(arrayData, extent,path):
    """
    write the given array data to the file 'path' with the given extent.
    """
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    metadata = driver.GetMetadata()
    if metadata.has_key(gdal.DCAP_CREATE) \
           and metadata[gdal.DCAP_CREATE] == 'YES':
        pass
    else:
        print 'Driver %s does not support Create() method.' % format
        return False

    # get rows and columns. or is it columns and rows?
    dims=arrayData.shape

    # could possible use CreateCopy from one of the input rasters...
    
    dst_ds = driver.Create(path, dims[1],dims[0], 1, gdal.GDT_Float32 )

    dst_ds.SetGeoTransform( [
        extent[0] , (extent[2]-extent[0])/dims[1], 0,
        extent[3], 0, (extent[1]-extent[3])/dims[0] ] )

    dst_ds.GetRasterBand(1).WriteArray(arrayData)

    return True
